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Abstract. The Procrustes distance is used to quantify the similarity or dissimilarity of (3-dimensional) 
shapes, and extensively used in biological morphometries. Typically each (normalized) shape is represented 
by N landmark points, chosen to be homologous (i.e. corresponding to each other), as much as possible, 
and the Procrustes distance is then computed as inf^j yZj—i \\Rxj — a;^- 1 1 2 , where the minimization is over 
all Euclidean transformations, and the correspondences xj <-> x'j arc picked in an optimal way. 
This (discrete) Procrustes distance is easy to compute but has drawbacks — representeding a shape by 
only a finite number of points, which may fail to capture all the geometric aspects of interest; a need 
has been expressed for alternatives that are still computationally tractable. We propose in this paper the 
concept of continuous Procrustes distance, and prove that it provides a true metric for two-dimensional 
surfaces embedded in three dimensions. The continuous Procrustes distance leads to a hard optimization 
problem over the group of area-preserving diffeomorphisms. One of the core observations of our paper 
is that for small continuous Procrustes distances, the global optimum of the Procrustes distance can be 
uniformly approximated by a conformal map. This observation leads to an efficient algorithm to calculate 
approximations to this new distance. 



1. Introduction 



Procrustes distances are used to compare shapes and quantify their (dis)similarity. In several applications, 
such as geometric morphometries [13], the shapes to be compared are continuous surfaces, on each of which 
homologous landmark points are selected, equal in number. The dissimilarity or distance between the surfaces 
iS and S' is then computed as the Procrustes distance between their corresponding landmark sequences 
X = (xg)p =1 and X' — (x' e )f =1 , which is defined as follows. 



Definition 1.1. Given two finite sequences X = (xi)™ =1 , X' — (a^)" =1 in R d , of equal length, with centroids 
— x' , and centroid sizes Sx, Sx'i respective!^] the classical Procrustes distance dp(X, X') between X and 



x 



X' is defined by 



inf 

Ren 



Rxi 
~Sx~ 



3_ 
Sx> 



1/2 



(1.1) d P {X.X') 
where 1Z is the group of Euclidean transformations (reflections, rotations, and translations) 



In some applications, it may be useful to consider weighted Procrustes distances, in which each label i € 
{1, . . . , L} can be given its own weight Wi in the computation of the centroids, the centroid sizes and the 
distance dp(X, X'); such weighting can be used to compensate, if desired, for possible imbalances in the 
distribution of the landmark points, when they occur more densely in some areas than in others. We shall 
assume in what follows that no such adjustment is needed, i.e., that the landmark points are considered 
(more or less) uniformly distributed. The normalization by the centroid size allows comparison of shapes 
irrespective of their scale. To achieve this without a normalization step, one would need to extend 1Z to 
the larger group of similarities, incorporating the (uniform) dilations as well. Note that other geometric 
extensive quantities could be used to normalize, with a very similar effect. 



The centroid of X is given by x = n 1 Xi ; the centroid size by Sx = [n 1 { x i ~ x ) 2 ] 1 
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The point sets X and X' are said to have the same shape if one can be obtained from the other by an 
appropriate combination of scaling, translating, rotating and (possibly) mirroring, i.e., if there exists R G 1Z 
and a G K such that aRX — X' . ("Shapes" of finite sets of points can thus be considered as orbits of these 
point sets under the action of similarity operations.) It is not hard to check that the Procrustes distance 
dp(-, •) is a metric on shapes, in the sense that it satisfies, for all finite point sequences X, X' and X" in R d 
of equal length, 

1) d P (X,X') =d P {X',X), 

2) dp(X, X') > 0, and dp(X,X') = implies that X and X' have the same shape, 

3) d P {X, X") < d P {X, X') + d P {X', X"). 

In all the above, the points Xi and x\ are ordered, i.e. corresponding points in X and X' have the same 
index. The correspondence between entries of X and X' can be "encoded" by a bijective correspondence 
map C : X — > X' that maps each Xj G X to its corresponding x'j G X' . In terms of this correspondence map 
the Procrustes distance dp(X,X') can be written as 

(1.2) d P (X,X') = d P (X,CX) = inf l\RX - CXI := inf [ V \\Rxi - Cxt 

Ren Ren \ ^— ' 

\»=1 

this recasts the minimization as a search for the map in 1Z that best approximates C, insofar as its action on 
X is concerned. 

Using a finite set X of landmark points X£ G <S, £ = 1, . . . , L as a proxy for the shape of a surface S, 
and taking the value of dp(X, X') to express the dissimilarity of the shapes of surfaces S and S' has some 
drawbacks, however. First, this approach compares only small discrete subsets of points sampled from the 
surfaces and therefore ignores "most" of their shapes. Second, and more importantly, it requires the user 
to carefully select corresponding landmark points on the two surfaces prior to calculating the Procrustes 
distance between the two landmark point sequences. This distance depends heavily on the exact choice of 
the landmark points. In geometric morphometries, one seeks to remove some of the arbitrariness of these 
choices by picking landmark points that are believed to be homologous, i.e., to truly correspond to each 
other, based on evolutionary arguments. This type of selection of landmark points requires considerable 
specialized expertise, and in some cases, even experts do not agree. In addition, morphologists interested 
in studying function of e.g. teeth are interested in moving away from landmark selection, and in using 
geometric information that encompasses more global features. 

This situation has motivated researchers to suggest alternative methods to compute distances or dissimi- 
larities between the shapes of surfaces. Even when these methods are based on continuous concepts, their 
numerical implementation requires some type of discretization, and thus often involves again discrete point 
sets X and X' (typically of larger cardinality than in landmark-based distances). The resulting distance 
can then still be written in the same form as the right hand side of ( |1.2[ ), with the important difference 
that C is no longer assumed to be given a priori. Instead, the map C is assumed to be determined by the 
full geometry of the surfaces S and S'; in practice, it has to be derived from the data themselves, meaning 
that both the correspondence C and the Euclidean transformation R must be determined numerically. (One 
could imagine a similar situation in the discrete case, if two sets X and X' were given, each with n points, 
without a correspondence map. In that case, a reasonable approach might be to select the map C : X — > X' 
for which \\CX - X'\\ is smallest.) 

A prominent method of this type is the Iterative Closest Point (ICP) algorithm [3]. This method alternates 
between determining C and R: the correspondence Ck ■ X — > X' is taken to associate to each point x' G X' 
the point(s) in X for which the image Rk-ix, under the best rigid alignment of X and X' obtained in 
the previous iteration, is closer to x' than to any other element in X'; the rigid alignment Rk is then the 
transformation R G 7Z that minimizes the distance \\RX — CkX\\. This algorithm is simple and robust, but 
suffers from several drawbacks. It may converge to a local minimum rather than the desired minimizer; 
this means that the limit may depend on the choice of the initial correspondence map C\ or the initial rigid 
alignment i? , whichever is picked to start off the algorithm. Of more concern is that the space of possible 
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correspondences C : S — > S' considered by the algorithm consists only of compositions of rigid motions 
and closest neighbor maps. This space of maps often contains high-distortion and discontinuous mappings, 
as illustrated in Figure [T] for the 1-dimensional situation; it also does not include a sufficiently rich set of 
diffcomorphisms (smooth bijective mappings). 




(a) The correspondence map C\ : X — > X' , in this case a length-preserving diffeomorphism. (S and S' each 
consist of 200 points, equispaced on the black horizontal line (X) and on the blue curve (X 1 ), respectively.) 




(b) Illustration of the Euclidean map R\ (moving the black line "up" ) and the correspondence map C2 ■ 

Figure 1. One-dimensional illustration of one step in the ICP algorithm. X 

and X' are point sets, each with 200 points, on two curves ("one-dimensional surfaces") 
<S (straight, black) and S' (wiggly, blue) of equal length, (a) The correspondence (in red) 
between X and X' that associates to each point x' € X' the point x <E X at the same 
arclength distance from the left end point of its curve, (b) Using this correspondence as 
an initial C\ : X — > X' , determine the Euclidean transformation (now a simple translation 
in the plane) R\ that minimizes \\RX — CiX|||, and move X to R\X\ the red lines now 
link each R\x £ R\X to the closest point in X' . The corresponding map C2 ■ X — > X' is 
discontinuous and highly distorting. 

Several authors have built extensions or generalizations of this approach, retaining the basic iterative princi- 
ple of ICP, interleaving the determination of correspondences Ck and transformations Rk in successive steps. 
Rangarajan et al. 15J formulate a variant on the Procrustes distance between two discrete sets of points in 
which the correspondence maps are unknown a priori. Their algorithm alternates between calculating opti- 
mal rotations and determining correspondence maps (bi- measures). For every fixed rotation i?, it computes 
the "measure coupling numbers" M,j from one point set to another, minimizing the average of the squared 
residuals j Mjj \\Rxi — x'j |, under the (soft) constraint that {Mij) i . =1 n is indeed a measure coupling. 
As is the case with ICP, this algorithm can still converge to a local rather than a global minimum, and the 
correspondence maps can still be "discontinuous and/or distorting". Ghosh et al. [5] use a similar framework 
(although not related to Procrustes or any other distance) with a smooth surface deformation mechanism 
together with closest point maps to determine both the correspondence maps and the transformations in an 
alternating iterative procedure. The algorithm in [8] requires user initialization (which may influence the 
outcome) ; the way correspondences are assigned can lead the deformation mechanism to ultimately produce 
a distorting and/or discontinuous map between the surfaces. 
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A common characteristic of the algorithms mentioned above, which often (in the limit or in intermediate 
stages) lead to discontinuous or distorting correspondence maps, is that the space they explore (implicitly 
or explicitly) to build correspondence maps is insufficiently rich in smooth bijections. 

In this paper we generalize the discrete Procrustes distance to continuous surfaces; in this formulation we 
use only smooth correspondence maps. Our construction leads to a non-linear functional over a huge and 
non-linear space of possible maps that we call the continuous Procrustes functional. Direct optimization, 
over its huge domain, of this generalized Procrustes functional is not feasible; we suggest that for many cases 
of interest, a different optimization suffices, over a (relatively) much smaller (and managable) subset of all 
possible maps, consisting of conformal mappings combined with specific area-preserving maps. One of the 
main results of our paper is a proof that the class of conformal maps uniformly approximates the globally 
optimal correspondence map in the regime where the continuous Procrustes functional takes on small values. 
Note that our approach thus provides a glimpse of the global minimizer (for the case of small continuous 
Procrustes distance) of a functional for which it is not known, in general, how to approximate the global 
minimizer in polynomial time. 

In addition to the theoretical constructions, we also provide an efficient algorithm, without user interaction, 
to construct (an approximation to) the continuous Procrustes distance and the optimal correspondence 
map between two surfaces, again in the case where this distance is small, i.e. where the surfaces are not 
too dissimilar. In practice, this algorithm performs very well, and is sufficiently fast to be used for the 
computation of pairwise distances for all pairs in reasonably large collections of surfaces (~ 100); see [TO], a 
first presentation of the main results of this paper at a workshop in June 2010, as well as [3], which uses the 
algorithm explained here in detail for three biological data-sets. 

A similar combination (conformal mappings composed with area-preserving maps) is used in the recent 
paper by Dominitz and Tannenbaum |5 , to construct good mappings from surfaces to a Euclidean spherical 
domain. The goal of [5 j is different, however; rather than seeking to define a distance between surfaces, 
that can be used for shape alignment, [5] is concerned with building a low distortion map from a surface to 
Euclidean domain, the inverse of which can then be used as a good parameterization for the surface. 

The paper is organized as follows. In section 2, we introduce our definition of the continuous Procrustes 
distance for homeomorphic 2-dimensional compact surfaces embedded in R 3 ; it involves a minimization that 
is unfeasible in practice. In section 3, we show that we can construct approximations of this distance by 
minimizing over appropriate perturbations of conformal mappings, which is much more tractable. In section 
4, we give the corresponding numerical algorithm and illustrate them with a concrete example. 



2. The Surfaces Procrustes Distance. 

Consider two homeomorphic compact 2-dimensional surfaces S, S' embedded in R 3 , endowed with the stan- 
dard metric induced from their embedding. Because the biological applications that motivated this work 
require comparing shapes irrespective of scale (see [13] and reference therein), we are interested in defining 
a scale-independent distance. We shall therefore assume that the two surfaces are normalized to have unit 
volume (area): 



dvoLj = 1=1 dvo\s> . 

IS JS' 

In building "good" correspondence maps between the continuous surfaces, we will be guided by what repre- 
sents a "good" correspondence between discrete (relatively small) sets of points that "represent" the surfaces 
when standard Procrustes distances are used. 

As mentioned in the introduction, great care is typically taken in the choice of sample points on surfaces that 
will then be used to compare the shapes of these surfaces. Landmark points on e.g. teeth or other bones are 
chosen so that they are homologous, i.e. "equivalent" from an evolutionary point of view. We are aiming 
for a landmark-free method; information of this type will thus not be available. Instead, we can use only 



geometric information given by the surface itself. Note that in ( 1.1 ), the different points Xi all play an equal 



role. When choosing discrete sets X, X', each consisting of n points, on the surfaces S and S' to represent 
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their respective "shapes", with the purpose of using them in a Procrustes distance calculation (1.1 1, it seems 
therefore reasonable to pick the points so that each represents an equal "share" of the surfaces; we shall 
interpret this here as representing an equal portion of the area of the surfaces. A correspondence map C 
that maps each Xi to its partner x\ can thus be interpreted as mapping portions of area 1/n of S to the 



corresponding portions of S' that have equal area; on the other hand, the sum in ( 1.1 ) can be viewed (up to 



a normalization) as a Riemann sum approximation to the integral of \\Rx — Cx\\ over S. 

This analysis suggests the following "continuous analogue" of the discrete construction. To involve the whole 
surface S (instead of just a set of sample points), we take C to be an area-preserving map from S to S'; for 
each fixed area-preserving C, we then define 

(2.1) d P (S,S';C)= M^(J \\Rx-Cx\\ 2 dvo\ s (x)^ ' . 

In the absence of landmark-type or other user-guided information we have to select C based solely on 
geometric information. Taking our cue from the discrete case, we want, given points (xi) i=1 on 5, to 
determine {x' i ) i=1 n on S' so that each x\ corresponds "as well as possible" to Xi. In other words, this 
suggests that C be picked so that S' and CS are optimally aligned, and that the continuous Procrustes 
distance be given by the corresponding value of dp(S,S';C). More explicitly, defining A(S,S') to be the set 
of diffeomorphisms (smooth bijective maps with a smooth inverse) from S to S' that are area-preserving, we 
set 

1/2" 



(2.2) n cP (s,s') 



inf 


inf { 


Js 


C£A(S,S') 


Ren \ 





\Rx — Cx\\ 2 dvo\s(x) 



In the remainder of this section, we establish several properties for the quantities defined in (2.1) and (2.2 1, 
establishing, among other results, that D c p(-, ■) defines a metric distance. 



We start by proving that the minimum in (2.1 1 is always attained 



Proposition 2.1. Given two homeomorphic surfaces S, S 1 of unit area, and an area-preserving map C from 
S to S' , there exists a rigid motion R* 6 1Z minimizing the functional J s \\Rx — Cx\\ 2 dvo\g(x) . 

Proof. Let R n £ TZ be a sequence such that 



J \\R n x-Cx\\ 2 dvo\ s (x)j -d P (S,S';C) 2 < ~. 



Let us represent each rigid motion as a composition of an orthogonal transformation and a translation: 

R n x U n x -\- t n , 

where U n £ M 3x3 and t n £ M 3 . Thinking of R n = (U n , t n ) as a vector in K 12 it is clear that there exists some 
compact set A C M 12 such that R n £ A for all n. Indeed, the orthogonal group 0(3) in its representation 
as a 3 x 3 matrix group is a compact set, and for sufficiently large n, the t n will all lie within some ball, i.e. 
||i n || < M for some M. Hence there exists some rigid transformation R* = (U* , t* ) such that, up to extracting 
a subsequence, R n — > R* as n — > oo. Lastly, R* realizes the infimum dp(S,S';C) since || 
\\U*x + t* — Cx\\ for every x £ S, and similar arguments as above imply that \\U n x + t n — Cx\\ is bounded 
uniformly in n and x £ S. The result then follows from the dominated convergence theorem. □ 



The following two propositions provide closed form solutions for (2.1|; their proofs follow the discrete case 
[6j in a rather straightforward manner. Note that we use that C is area-preserving to establish these formulas 
(for the translation part). First, we show that the translational part t* takes the centroid of S to the centroid 
of 5': 

Proposition 2.2. IfS andS' both have their centroids at the origin, i.e. J s xdvols(x) = = J s , y dvol^i (y), 
then the translational part t* of the optimal rigid motion R* is zero. 
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Proof. Assume the surfaces S, S' are centered as described in the assumptions of the theorem. Differentiating 
J s \\Ux + t — Cx\\ 2 dvo\s(x) with respect to each of the coordinates of the vector t and plugging in U = U* 
and t = t* , we get 

= 2^ (U*x + t* -Cx)dvols(x). 

Rearranging the above equality and remembering that S and S' have unit area, we get 

t* = [ Cx dvo\ s (x) - [ U*x dvo\ s (x) 
Js Js 

= [ y dvol s >(y) - U* I x dvo\ s (x) = 0. 
Js 1 Js 



□ 



Next, the orthogonal transformation part: 

Proposition 2.3. If S andS' both have their centroids at the origin, i.e. J s xdvols (x) = = f s , ydvols>(y), 
then the optimal orthogonal transformation U* can be written as 

U* = WQ T , 

where W,Q are the orthogonal transformations from the Singular Value Decomposition (SVD) 



Is 



x{Cxf dvol s {x) = QSW 1 



where S = diag (<ti, 02, 03) is a diagonal matrix with the singular values of U* on the diagonal. 
Proof. Expanding J s \\Ux — Cx\\ 2 dvo\s(x) we get: 



/ ||C/a;-Ca;|| 2 dvol5(a;) 
Js 



\\x\\ 2 dvo\ s (x) - 2 / x T U T Cx dvo\ s {x) + / \\Cx\\ 2 dvol s {x), 
is Js J 

where we used that U T U — Id. The sought-for U* therefore must maximize 



Note that 
and therefore 



E{U) = J x T U T Cx dvols(x). 
x T U T Cx = tr (x T U T Cx) = tr (llx (Cx) T ) 



E(U) = tr 



U J x(Cxf dvo\ s (x) 



tr [W T UQS] = tr 



US 



tr [UQSW T ] 



where U := W T UQ, and we used the SVD decomposition. Note that the last term cannot be greater than 
o~i + 02 + 03 since all entries of the orthogonal matrix U have absolute value at most 1. Further note that 
taking U = WQ T achieves this upper bound. The uniqueness is also clear. □ 

We now prove: 

Proposition 2.4. For each fixed area-preserving map C from S to S' , we have 

(1) d P (S,S';C) >0 

(2) d P (S,S';C) = dpiS'&C- 1 ) 

(3) dp(S,S';C) = implies that S and S' are congruent. 



THE CONTINUOUS PROCRUSTES DISTANCE BETWEEN TWO SURFACES 



7 



Moreover, if S" is a third surface, and C is an area-preserving map from S' to S" , then 
(4) d P (5,5";C'oC) < d P (S,S';C) + d P (S',S";C). 



Proof. First, it is clear that dp(S,S';C) > 0. If dp(S,S';C) — 0, then we know by Proposition 2.1 that 
there exists a rigid transformation R* such that 



J \\R*x-Cxf dvol s (x) = 0. 



Since we are dealing with smooth surfaces S, this implies that \\R*x — Cx\\ = for all x € S; since the range 
of C is all of S' (because C is a bijective diffeomorphism) it follows that R*S and S' are equal as sets, so S 
and S' are congruent. 

Next we prove symmetry: 

dp(S,S';C) 2 = inf f \\Rx - Cx\\ 2 dvo\ s (x) = inf f {{RC^y - yf dvo\ s ,{y) 
J s Ren J s> 

= inf f \\C- 1 y-Ry\\ 2 d V ol s ,(y)=d P (S / ,S;C- 1 ) 2 , 

Ke/,i JS' 

where the second equality uses the fact that C is area-preserving. 
For arbitrary R G 1Z, we have 

/ r \ V 2 

d P {S,S";C oC) = inf / \\Rx - C o Cxf dvo\ s (x) 



< inf 
Ren 



= inf ■ 

Ren 









Rx — RCx 


W,) V2 } + (/ s 


RCx — C' o Cx 









Rx — Cx | dvolg {x) 



= d P (S,S';C) 



Ry - C'y 



1/2* 



dvols* (y) 



Ry - C'y 



dvolg (x) 

1/2 



1/2 



d-vo\ S ' (y) 



1/2 



By taking the infimum over all R € 1Z we obtain the desired result. 



□ 



Having established these properties for dp(S,S';C), we can now minimize this over C € A(S,S'), i.e. we 
have 



(2.3) 



T> cP (S,S')= inf d P (S,S';C); 



if the infimum is achieved by some C* S C, we declare this C* to be our desired correspondence map. 

Whether such a minimizer exists is a delicate question, and we do not have a proof or counter example 
for the general case. However, if we restrict the class of maps C to bi-Lipschitz maps with some a priori 
bound on the maximal dilation, then such a minimizer does indeed exist; moreover this minimizer is also 
bi-Lipschitz with the same bound. 

Proposition 2.5. For arbitrary B > 0, let Bb{S,S') be the set of bi-Lipschitz area-preserving diffeomor- 
phisms from S to S 1 such that, for all x, y G S , B^ 1 &s{x, y) < dgi(Cx,Cy) < Bdg(x,y), and let d$, 
resp. d$i denote the geodesic distances on S, resp. S' . Then there exists a minimizer in Bb(S,S') for the 
restriction to Ss(5,5') of the functional D c p(S,S'; •). 



Proof. It is straightforward that Bp(S, S') is a closed subset of C(S, S'), the set of continuous functions from 
S to S' , equipped with the topology of uniform convergence with respect to ds and ds'. By the definition 
of Bb(S,S'), the functions in Bp are equicontinuous. It then follows from the Ascoli-Arzela theorem for the 
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continuous functions on compact metric spaces that Bb(S,S') is compact. 

It is also easy to see that the functional dp (5, 5'; •) is continuous with respect to the topology of uniform 
convergence on C(5,5'). It follows that the restriction of dp(5,5';-) to Bb{S,S') is a continuous map 
from a compact space to R. Let now (C„) ngN be a minimizing sequence in Z3b(5,5'), i.e. dp(5,5';C„) — > 
ml CeB B (5,S') dp (5, 5';C) as n — > oo. By the compactness of 23^(5,5'), the sequence (C„) ngN has a uniformly 
converging subsequence ; if we denote its limit by C * , then it follows that C * E Bb (S,S'), and dp(S,S';C*) = 
ini ceB B (s,s')dp{S,S';C). □ 

We note here that all the further proofs and results in the paper will remain valid (mutatis mutandum) if 
we replace everywhere the class of general area-preserving maps by the more restricted class of bi-Lipschitz 
area-preserving maps. 

Even when the existence of a minimizer is not guaranteed, it is possible to prove that D cP (5,5') defines a 
metric up-to congruence relation: 

Theorem 2.6. D c p(5,5') defines a metric between surfaces wp -to -congruence, that is, D c p(5,5') > 0, 
D cF (5,5') = D cF (5',5), D cP (5,5') < D cP (5,5") +T> eP (S",S'), and D cP (5, 5') = only if S and S' 
are congruent. 

Proof Clearly D cP (5,5') > 0. 



(i*n) 



If T) c p(S,S') = then we have a sequence (C„)„ s n C A(S,S') and (by Proposition 2.1| a sequence 
C 1Z of rigid motions such that: 



2 1 

\\R n x-C n x\\ a\o\ s (x) < 



s 



n 



3 ' 



By extracting a subsequence [Rlj (with R* k := R nk ), we can assume that \\R^ — R**\\ — > 0, with 
R** e TZ as k -> oo, and that (with C| := C nj J 



y ||i?**.T-C^|| 2 dvol 5 (x) < p. 



Set now B nJ := {x ; - C*ir|| 2 > 1/ff, 4 £ := {x ; lim sup^^ - C£a;|| 2 > 1/ff = n r „ eN U„> m 

S„,£ , and v4oo := Ut>iAg — {x; limsupj.^^ — C%z\\ ^ 0}. Since Ag C for all i, it follows 

that vol 5 (^l 00 ) = lim^oo vols(j4^) . We have vols(.B n ^)< I J s \\R**x - C^x\\ 2 dvo\ s (x) < i/n 3 , hence 
vol 5 ( U„> m < i/m 2 and thus vol s (A^) = vol 5 ( n meN U n > m B n ^ < M meN £/m 2 = . It follows that 



vol s (A 00 ) = 0, or vol s ({a;; limsup^^ \\R** x - C* k x\\ ^ 0} 



0. 



Therefore R**x — lirrife-^oo Ck{x) for x e 5 \ A^, implying R**(S\ A^) C 5' = 5'. Since every open disk in 
5 (with respect to the geodesic distance on 5) has area strictly greater than in 5, 5 \ Aoo is dense in 5. 
By the continuity of R** it follows that R**(S) C 5'. 

Let's assume now (hoping to derive a contradiction) that there exists a point y £ 5' such that y ^ R**(S). 
Since R**(S) is a closed set there must then exist a set O C 5' with positive area such that Or\R**(S) = 0. 
This is a contradiction since R** : 5 — > 5' is an isomctry and in particular area-preserving. Hence R**(S) = 
5', showing that 5 and 5' are congruent. 

Symmetry is easy to establish as follows: 

D cP (5,5') = inf dp(S,S';C)= inf d P (5', 5; C _1 ) 

= inf dp(5',5;C)=D cP (5',5), 
ee-A(S',s) 

where we used that C e -4(5,5') iff C" 1 e A(S',S). 
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Lastly, for the triangle inequality, we have, by Theorem |2.4| for every C G A(S, S') and every C G A(S',S"), 
B cP (S,S") = inf d P (5,5";C") < d P (S,S";C o C) < d P (S,S';C) + d P (S',S";C). 

C"GA(S,S") 

Taking the infimum over C G A(S, S') and C G A(S' , S") we get 

D cP (5,5")< inf dp(S,S , ;C)+ inf d P (S',S";C) = T) cP (S,S') +T) cP (S',S"). 

□ 

We conclude this section by providing an approximation result: given two surfaces S, S' and a correspondence 
map C G A(S, S'), we would like to approximate the centroids J s x dvols (x), J s , y dvols' (y) and the integral 
defined in Proposition [23J these approximations will be used to compute approximations to the optimal rigid 
transformations and to the distances d P (S, 5'; C). To this end we will use a simple rectangle-type integration 
formula that we describe now. Let Q — {qe}g = i C S be a set of points such that their corresponding Voronoi 
cells {T|}^ have approximately equal surface area; in practice, such a set of points can be determined by 
means of the Furthest Point Algorithm (FPS) [7]. Using the notation vols (Ye) — J T ^ dvols we then have 



(2.4) 



/ f(x)dvol s (x) & V/fe)vol 5 (T,). 
Js 



The error made in this approximation can be estimated in terms of the fill distance t]s(Q) of the set Q, 
defined as 



(2.5) 



r/ s (Q) : = sup jr G K 3x G M s.t. B s (x,r) f)Q = $\, 



where Bs(x,r) — {q G S | ds(x,q) < r}, with ds(x, q) the geodesic distance on S between x and q. 
Intuitively, the fill distance f]s(Q) is the radius of the largest geodesic open ball that can be placed on the 
surface S without including any point of the (discrete) set Q. In other words it is the largest "circular hole" 
in the sampling Q C S. We have then 



Proposition 2.7. The error of the approximation (2.4) has the following upper bound. 



f(x)dvol s (x) - ^2f(q £ )vols(T e ) 



< sup \f(x)-f(q e )\<Mr)s(Q), 



where M is a bound on the norm of the gradient Wsf °f f ■ Hence the error is linear in the separation 
distance. 



Proof. Writing S = U^T^ we get 

f f(x)dvols(x)-Y J f(l^ols(T l ) <J2[ !/(*)-/(&) I dvol s (a) 

J S is f)>J If 



< sup \f(x)-f(q l )\Y^vol s (T l 
= sup \f(x) - f(q e )\ , 



where the last equality uses vol 5 (5) = 1. Now take arbitrary x G T e , and denote by j(t) : [0, d s (qe, x)] — > S 
the arc-length speed geodesic curve connecting qi and x. Then, 



rd s (qi,x) j rd s (qe,x) 

m f(m) = J o Jt if(i(t))] dt = (V 5 /(7(*)), i(t)) s dt- 
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Using the Cauchy-Schwarz inequality, 

i-d s {qi,x) 

\f(x) - f(q t )\ < sup \\V S f(x)\\ s / ||7(i)|| s dt 

xes Jo 

= sup||Vs/(x)|| s ds(qt,x). 
xes 

Lastly, the inequality d$(q£,x) < rj(S) can be derived directly from the properties of Voronoi cells (see 
Lemma D.2 in [H]). □ 

3. MOBIUS TRANSFORMATIONS AS A REDUCED SEARCH SPACE 

Computing the surface Procrustes distance, as we defined it above, amounts to solving a hard optimization 
problem: unfortunately, the sets A(S, S') or Bb{S,S') are formally infinite dimensional manifolds, and 
therefore extremely hard to search in practice. Our key idea is to replace the search space A(S,S') in the 
variational formulation ( |2.3| by another, much smaller, set of maps. The core observation is that the set of 
conformal (or anti-conformal) mappings between S and 5', which has a finite (and small) dimensionality, 
gets "close" (in some sense to be made precise below) to the minimizing C £ A(S,S'). In particular, we shall 
see that if T) c p(S,S') is small, then the minimizing area-preserving map C in (|2.3[) is close to conformal. 



Let us explain this in some more detail. We are particularly interested in computing (approximate) con- 
tinuous Procrustes distances for "close" pairs 0j. In those cases the insight that (close to) optimal C have 
to be close to conformal leads us to a strategy that involves minimizing over a much smaller set of maps. 
To achieve this, we shall make use of a nonlinear procedure Pr that "transforms" a map that is close to 
.4(5,5') into an area-preserving map, i.e. to an element of A(S, <S'). This nonlinear transformation leaves 
elements of A(S, S') unchanged, and can thus be interpreted as a nonlinear projection procedure (hence the 
notation). The smaller set of maps over which we shall minimize is then the image in A(S,S') of the family 
of conformal maps (from S to S'), transformed by Pr. 

As a search space, the family of conformal mappings is a much more "friendly" setting than A(S,S') or 
Sb(5,5'). First, by the uniformization theorem the conformal (or anti-conformal) bijective mappings can 
be characterized completely, and an explicit parameterization can be given in terms of a small number of 
parameters. For instance, the family of conformal bijective mappings between two disk-type surfaces S, S' 
is represented by (disk-preserving) Mobius transformations. Each mapping in this family is completely 
characterized by 3 real (bounded) parameters; therefore the search over the space of conformal mappings 
can be done efficiently. Second, Mobius transformations are smooth bijective diffcomorphisms, so that our 
candidate search space consists of only "nice" intrinsic mappings. 

To motivate why we would consider restricting ourselves to conformal mappings (or their deformations 



through Pr) for the optimization, we note that for S, S' such that T) c p(S,S') = 0, the infimum in (2.3) 
is achieved for some R S 1Z (by Theorem |2.4[ ); in this case the minimizing C = R is obviously conformal. 
However, we prove below the stronger result that a correspondence C : S — > S' for which the distance 
dp(S 1 S';C) is small can be approximated (under rather mild assumptions on the regularity of C) by a 
bijective globally conformal mapping from S to S' . 

We start with a few simple lemmas. The first Lemma is proved in |16j : 

Lemma 3.1. Let S C K 3 be a compact 2-manifold with the induced Riemannian metric g. Then 

d s (x,x')- \\x-x'\\ | < C s d s (x,x') 3 , 

where ds(x,x') denotes the geodesic distance between x and 'I! denotes the Euclidean distance 

between these points, and Cs depends only on the curvature ofS. 



Next, we prove a result concerning the approximation of the norm of the differential of a map: 
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Lemma 3.2. Let S, S' C R 3 be compact 2-manifolds with the induced Riemannian metrics g,h ( respectively) . 
Let F : S — > S' be a smooth map, and denote by DF X the differential of F at arbitrary x G S . Then for 
5 > sufficiently small, the following holds: for all x G S, there exists x' in the boundary dB g (x,S) of the 
S-radius geodesic ball centered at x, dB g (x,5) := {u G S\ ds(x,u) = 5}, such that 



\\F(x)-F(x>) 



\\x- 



\DF. r 



<g,h 



< C K 6, 



where \\-\\ gh is the operator norm associated with the norms ||-|| and \\-\\ h in the usual way, i.e. \\L\\ gh = 
su Pf^o • Here, C K depends only on the maximum of the surfaces' curvature and norms of second order 

differentials of the mapping F . 



Proof. For £ G R 2 , we set V(5) = {£, G R 2 | ||£|| < S}. Fix x and take 6 > small enough so that the 
exponential map exp x : T>{25) —> B g (x,25) is a difFeomorphism. For every x' G B g (x,S), we denote by the 
vector ^1 the vector in T>{5) such that exp x (£ x i) = x'\ in other words, £, x > is tangent to the geodesic on S 



that goes from x to x', and \\£ x 
j(t) = t£ x >, < t < 1. Then 



ds(x,x'). Denote F — F o exp x : T>{8) — > S', and consider the line 



F(x') - F(x) = 



F(j(t)) dt= / DF 7{t) j(t)dt 
J Jo 



dt 



(3.1) 



J (dF i(0) + Ofll&'ll))^'* = DF X & + 0(d s (x, x' f ) . 



Let i e dV(5) be such that DF X (£) 



DF X 



g,h 



|^|L. Remember that at exp 1 {x) the pull-back metric 



tensor (exp* x g) equals and therefore ||^|| 



Also note that DF X (£) 



DF X (£) 



since the 

h 

metric h is induced by the ambient Euclidean metric of R 3 . Now set x' = exp x (^), and take the Euclidean 

OF, d s (x,x') + 0{d s {x,x') 2 ), 



norm of both sides of (3.1 ). Then we have 

||F(^)-^)|| = 
and therefore 



g,h 



\\F(x')-F(x) 



d s (x,x') 

Using Lemma [3T] this leads to the desired estimate 



DF*t 



g,h 



0(d s (x,x')). 



a 



Next, we define the cone condition for a surface S: 

Definition 3.3. We say that a compact 2-manifold 5cR 3 satisfies the (a,6)-cone condition, where a > 
and 9 G (0, 2tt], if for every i£5 there is a unit vector n in the tangent plane T X S such that the exponential 
map exp x is well-defined on the cone cs(<J,0;n) = {£ G M 2 ; \\£\\ < a , (£, n) > ||^|| cos(0/2)}, and is one-to- 
one on the whole cone. 

In other words, the surface S satisfies the (a, #)-cone condition if for every i?5, there is a "fan" , spanning 
at least an angle 9, of geodesies that leave x and continue, within S, for at least a distance a (w.r.t. the 
metric induced on 5 by I 3 ), without intersecting themselves or any other geodesic in the fan. We have now 

Lemma 3.4. Let S C 1R 3 be a compact 2-manifold satisfying the (a,9)-cone condition. Then there exist 
constants p , T > depending on a , 9 and on the curvature k of S such that for all u G S and all r < p, the 
area of {x G S\ \\u — a:|| < r} is bounded below by Tr 2 (with ||-|| standing for the Euclidean norm in M?). 
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Proof. By Lemma 3.1 there exists a constant R > (depending only on the curvature of S) such that for 
all x, y £ S satisfying ds(x, y) < R we have 

ds(x,y) > ^ \\x-y\\ • 

Set po = 2 min {R, a}, and fix an arbitrary u £ S. We have then, for all r < p , that {x £ S\ ds(u, x) < r/2} C 
{x £ S\ \\u — x\\ < r}, and consequently 



dvols(x) > / dvols(x). 

{x£S\ ||u-x||<r'} J{x£S\ ds(u,x)<r/2\ 

Now introduce polar coordinates (r, (f>) on the tangent plane T U S, so that the vector n (with respect to which 
the cone condition holds at u) is aligned with the direction = 0. With respect to this coordinate system, 
the exponential exp u maps [0, a] x [—8/2,9/2] to S, and the metric density can be written as |17j : 

V^(r,0)=r-r 3 '^ ) + O (r 3 ). 
Since r < p < 2a, the sector [0, r/2] x [—0/2,9/2] is contained in [0, a] x [—9/2,9/2] and we have 

9/2 ,r/2 

dvol 5 (a:) > / / y/g(r,(j))dTd4> 



{xeS\ ds(i,j)<r/2} J -0/2 JO 



8 4 k(u)9 



2 U 4 

r r 



+ o(r 4 ) =r 2 l(l+0(r>)), 



384 

where, as usual, the absolute value of the 0(r 2 ) term is bounded above by Cr 2 , for some C > 0, for all r 
smaller than some r\. Setting p = min( ( oo, fi, l/v / 2C) and T = [1 — min(l/2, C rf)] 9/8 we obtain, for r < p, 

dvo\ s {x) >Tr 2 , 

{xeS\d s (x,y)<r/2} 

completing the proof. □ 



We are ready to prove the main result of this section, which provides a bound on the conformal distortion 
disc of the optimal area-preserving "alignments" C for surfaces S and S' that are close to each other in the 
continuous Procrustes distance. The conformal distortion disc(x) of C at x is defined as the ratio between 
the two singular values of the matrix obtained by expressing the differential DC X with respect to orthonormal 
bases in T X S and Tq x S' , respectively. 

Theorem 3.5. Let S, S 1 C K 3 be 2-manifolds with induced Riemannian metrics g, h (respectively), with 
curvatures bounded above by k, and satisfying the (a,9)-cone condition. We consider area-preserving diffeo- 
morphisms C : S — > S' with first and second order differentials bounded by M . Then, for sufficiently small 
e, the bound dp(S,S';C) < e implies the following bound on the conformal distortion disc of the map C: 

sup dis c (x) < 1 + 0(e 1/4 ), 

x£S 

where the constant in the O-notation depends on only k and M . 



Proof. Denote by R £ 1Z the rigid motion for which the infimum in (2.1 ) is attained for C 



The first step in our proof is to derive a uniform bound on ||-R(x) — C(x)||. We start by noting that the 
function q(x) = \\R{x) — C(x)\\ is Lipschitz with a constant A dependent only on M. Indeed, we have 



\\R(x)-C(x)\\-\\R(y)-C(y)\[ 



< \\R(x)-R(y)\\ + \\C(x)-C(y)\\ 
= \\x-y\\ + \\C(x)-C(y)\[ . 



By assumption, ds>{C(x),C(y)) < Mds(x,y). By Lemma 3.1, dg(x,y) < 3/2 — y\\ if d$(x,y) is suffi- 
ciently small. On the other hand, we have, for all x' , y' £ S' , \\x' — y'\\ < dg> (x' , y'), dgi is the metric 
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induced on S' by the Euclidean metric in M 3 . Thus \\C(x) — C(y)\\ < 3M/2 \\x — y\\ when ||x — y\\ is suffi- 
ciently small. Since on the other hand S' is compact and thus bounded, \\C(x) — C(y)\\ is bounded uniformly 
in x, y, regardless of ||a; — y\\. ft follows that there exists a constant A, depending only on the geometric 
properties of the surfaces S and S' , and on M, such that 



\R(x)-C(x)\\-\\R(y)-C(y)\\ 



<Hx-y\\ 



Suppose q attains its maximum Q in it G S. Set a — min (p, , with p > as in Lemma 
must have 



3.4 



Then we 



max(0,Q- A||u-x||) 2 dvol s {x) < / \\R(x) - C(x)\\ 2 

{x£$\ \\u—x\\<a} J {xGS\ \\u—x\\<a} 

< d P (5,5';C) 2 < e 2 . 

On the other hand, we also have, by Lemma |3.4| and using Q — X \\x — u\\ > Q — \a on {x G <S| |ju — x|| < a}, 

max(0,Q- A||u-a;||) 2 dvo\ s {x) > / (Q - Xa) 2 dvo\ s {x) 

{x^Sl \\u—x\\<a} J {x£$\ \\u—x\\<a} 

2A 



(Q - Xa) 2 I dvol s (x) > ^-Ta 2 = — Q 2 min ( p 

{xeS\\\u-x\\<a} 



This implies, in particular, that 

r 2 / g\ 2 , 2 

— Q mm ( p, — ] < e 



2A/ 

If Q/(2X) > p, then it follows that TQ 2 p 2 /4 < e 2 , hence (by using Q/(2A) > p once again) TA 2 p 4 < e 2 . 
Note that T, A and p are constants that depend on only the geometrical bounds that we impose on S, <S' 
separately; a priori they bear no relationship to whether or not the continuous Procrustes distance between 
the surfaces is small. With a left hand side independent of e and strictly positive, the inequality above can 
therefore not be satisfied if e is sufficiently small; more precisely, if e < T 1 / 2 A p 2 ), then this case is excluded. 

For sufficiently small e, we have thus Q/(2A) < p, implying rQ 4 /(16A 2 ) < e 2 , or Q < 2 A 1 / 2 T" 1 / 4 e 1 / 2 . In 
other words, there exists a constant C\ > 0, dependent on only A, p, M and k, such that, for sufficiently 
small e, 

max||JZ(a:)-<7(a:)|| = Q < C± e 1/2 , 
which is the desired uniform bound on — C(a;)||. 



Second, by Lemma 3.2 we can take y G dB g (x, e 1 / 4 ) such that 

\\C(y)-C(x)\\ 



\\x-y\\ " x " g ' h 



DC x \\ a ,+0(e^). 



Using the triangle inequality as well as ||-R(a;) — R(y)\\ = \\x — y\\, and applying Lemma 3.1 we obtain 
\\C(y)-C(x)\\ \\C(y)-R(y)\\ + \\x-y\\ + \\R(x)-C(x)\\ 1/4 
\\x-y\\ ~ \\x-v\\ -l + ^e ), 

and thus 

pXy flj/ ,<l + 0(eV4). 

Lastly, since ||DCa;|| g h equals the larger singular value of the matrix for DC X w.r.t. orthonormal bases of 
T X S and T C ^S' (respectively), and since C is area-preserving (implying that the determinant of this 2x2 
matrix equals 1) the conformal distortion of C at x is H-DC^H 2 h , and thus 

disc(x) < l + 0(e 1/4 ). 

□ 
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Theorem 3.5 tells us that area-preserving diffeomorphisms associated to small surface Procrustes distances 
have small conformal distortion everywhere. We will next use the theory of quasi-conformal (QC) maps to 
see that, for disk-type surfaces, this implies that such maps then must be "close" to conformal maps. 

For the sake of convenience, we restrict our discussion here to the case of disk-type surfaces here (similar 
results can be shown for sphere-type surfaces). More precisely, we start with two disk type surfaces S,S' C 
K 3 with induced metric tensors <?, h (respectively), and we consider a global conformal parametrization 
(uniformization) of each onto their canonical domain, : S — > V, : S' — > T>. The surfaces are then 
intrinsically represented by their conformal factors fi and v. In other words, the push-forward metric tensors 
of S, S' under the maps Sfr' are given by (\E , *(7)[z] = /j,(z)dzdz, and (\f^/i)[i<;] = v(w)dwdw, respectively. 
The conformal factors also act as "density functions" in the sense that the area in S of an arbitrary Borel 
set C S can be written as vols(f2) = J^q-, fi(z)dxdy, where z = x + iy; similarly for the surface S'. 

Now every conformal mapping from S to S' can be written as vp'^ 1 omof, where m ranges over the Mobius 
transformations of the unit disk that preserve its boundary: 

■ w z-a 



(3.2) m(z)=e ] - 

where 9 € [0, 27r), a <E T>. This family of transformations has three degrees of freedom (one for the angle and 
two for the complex number a); we denote the family by Mob{T>). 

Likewise an area preserving (and orientation preserving) map C from S to S' can be "transported" to T> by 
means of 'J and leading us to consider instead C tI := oC o mapping T> to itself. We will use 

QC theory to show that, if dp(S,S';C) is small, then C tr is close to an element of Mob(D), with respect to 
the maximum norm over the unit disk T> — {z £ C \ \z\ < 1}, at least if C is orientation preserving. If it 
is orientation reversing, it is close to an anti-conformal map. We provide details below for the orientation 
preserving case; the reversing case is entirely similar. 

By an appropriate Mobius change of coordinates m, replacing $ by = m o <3>', we can even ensure that 
C tr := m o C tt has and 1 as fixed points. Abusing notation, and denoting C tr by C again, we thus assume 
C(0) = 0, and C(l) = 1. We shall show that C is close to the identity, which means that C tT is close to m _1 , 
and thus that the original area preserving map from S to <S' is close to the conformal map from S to S' 
given by o rhr 1 o 

We consider, as is very customary in complex analysis, derivatives with respect to z and z of the diffcrcntiablc 
map C from the subset I? of C to itself, i.e. 

dC dC dC , dC dC dC 

— = i — and — = — + i — , where z = x + iy. 

oz ox oy az ox oy 

We define the complex dilation of C by 

_ dC I dC 
® dz' dz 

For orientation preserving C we have (see for example, pQ): 

1 + \g\ disc - 1 

dlSc = THei ' = rf^TT' 

Theorem |3.5| therefore implies, uniformly on T>, 

(3.3) Q = 0{e 1 ' i ). 

We will use the following existence and uniqueness theorem for the Beltrami equation (see [9], Theorem 
4.30): 

Theorem 3.6. For every g : C — > C measurable such that WqW^ < I, there exists a homeomorphism f of C 
onto C which is a quasiconformal mapping ofC with complex dilation g. Moreover, f is uniquely determined 
by the following normalization conditions: /(0) = 0, /(l) = 1, and /(oo) = oo. 
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As is customary, we will call normalized solution any solution of a Beltrami equation that satisfies the 
normalization conditions. 



Theorem |3.6| requires the complex dilation p to be defined on all of C. Before applying it, we thus need to 
first obtain p on all of C, which we do by extending C from T> to the entire complex plane C by reflection: 



C(z) 



C(z) \z\<l 
l/c(l/z) |*| >1 



Note that the extension C is normalized, that is, it satisfies C(0) = 0, C(l) = 1, C(oo) = oo. Moreover, this 
extension preserves the conformal distortion, that is, for \z\ > 1: 

disp(z) = disc(l/z). 



It follows that this extension of disc to all of C still satisfies (3.3). We now have 

Lemma 3.7. The extension C : C — > C is the unique normalized solution to the following Beltrami equation: 



where g is a complex dilation (a.k.a. Beltrami coefficient) defined by 

g(z) = 



g(z) \z\ < 1 

\z\ = l 



0(1/2) (§*) ze\z\>l 

Proof. A straightforward calculation shows that C has the Beltrami co effici ent g almost everywhere. As 



mentioned above, C satisfies the normalization conditions of Theorem 3.6 and therefore the uniqueness 



follows from Theorem 13.61 □ 



We will next use Proposition 4.36 from [9], the statement of which is: 

Theorem 3.8. If WgW^ —> 0, then the normalized solution of the Beltrami equation f e converges to the 
identity in the maximum norm on V, \\f g — IdW^ —> 0, where Id(z) = z. 

The proof of Proposition 4.36 in [3] actually demonstrates a slightly stronger claim: 
Theorem 3.9. If 1 1 ,o| 1 00 — > 0, then the normalized solution of the Beltrami equation f Q satisfies 

||/ e -«|loc<^IML> 

on T>, for some constant M > independent of sufficiently small g. 
Combining Theorems |3.5| and |3.9| finally yields 

Theorem 3.10. Let 5, S' C K 3 be 2-manifolds with induced Riemannian metrics g, h (respectively), with 
curvatures bounded above by k, and satisfying the (a, 9) -cone condition. We consider area-preserving and 
orientation-preserving diffeomorphisms C : S — > S' with first and second order differentials bounded by M . 
Let \1/ : S — > T>. : S' — > T> be uniformizing maps of S, <S' onto the disk. Let m be a disk-preserving 
Mobius transformation such that / ^moip'oCo "J -1 : D — > D satisfies /(0) = 0, f(l) = 1. Then the bound 
dp(S,S';C) < e implies the following bound: 

ll/-W|loo = 0(6 1/4 ), 

where Id(z) — z is the identity map, and where the constant in the O-notation depends on only k, a, 9 and 
S. 
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Figure 2. A surface (left) and its density function over the uniformization disk. The local 
extremas of the conformal density function are shown as black dots. 

The orientation reversing C : S — > S' are close to the anti-Mobius transformations that can be calculated 
from the Mobius transformations by setting 

(3.4) m(z) = m{z), 

where m is any Mobius transformation. 

As described earlier, we use this theorem as a guide to build an efficient search algorithm to compute (an 
approximation to) dp(S,S') for surfaces <S, S' that are not hugely dissimilar. Since area-preserving maps 
C from S to S' that are close to minimizing dp (S,S';C) must be close to conformal, we start by searching 
Mob(S,S') to find the conformal or anti-conformal map m that minimizes dp(S,S / ;m). We then transform 
this m into a nearby area-preserving diffeomorphism by means of a nonlinear transform Pr, still to be defined 
below. We expect (but do no prove) that Prim) is then a good approximation to (nearly) minimizing C. 
Note that there are no guarantees that this approximation process, in which we replace A(S, S') by the proxy 
Pr (Mob(S,S')), preserves the triangle inequality property of D c p(<S,<S'); the approximations we compute 
therefore result in a measure of dissimilarity rather then a distance. 



4. Searching appropriate Mobius candidates and massaging them into area-preservation 



In the previous section we showed that it is useful to first find a Mobius transformation m for which 
dp(S,S';m) is small; we show in subsection 4.1 below a practical strategy for obtaining such candidate 
Mobius transformations that is fast and efficient for our applications. To obtain a better approximation of 



the optimal element of A(S,S') from these candidate m € Mob{S ,S'), we will, in subsections 4.2 and 4.3 



deform each of them into a nearby area-preserving map. That is, for every m £ Mob(S, S'), we will construct 
a map f m : S — > S 1 such that f m o m 6 A(S,S'). 



4.1. Searching the Mobius group. By Theorem |3.10| we know that an area-preserving diffeomorphism 
C : S —> S' that produces a small continuous Procrustes distance dp(S,S';C), is close to a Mobius trans- 
formation, when written in uniformizing coordinates. Hence, we first describe how we search for candidate 
Mobius transformations m € Mob(S, S') that (we hope) are already close to area-preserving for our appli- 
cations. As mentioned above, the Mobius group between two disk-type surfaces has three real degrees of 
freedom: prescribing the image wq G T> of one point zq e T>, as well as one angle 9 G [0, 2ir), uniquely defines 
a disk-preserving Mobius transformation m : T> — > T>. To speed up the search, we start by determining a 
mapping for which the density peaks, i.e. the local extrema of the density v(m(z)) (more or less) corre- 
spond to those of n(z). To that end we first extract, for each surface, a set of extremal points Is, Is' (local 
maxima and minima) defined by local extrema of the corresponding density functions fi, v, respectively. See 
Figure pi where the black points show these extremal sets. In practice, we find that, in the application (to 
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bone surfaces) that first motivated us, these points (intimately related to extrema of Gauss curvature) were 
likely to contain at least one pair of corresponding points, across a wide range of examples; this feature has 
presisted for other families of examples we examined. Note that this definition of Is,Is' is not invariant 
to Mobius transformations in the sense that for any Mobius transformation the extrema of the pulled-back 
density v(m(z)) \m'(z)\ 2 are not, in general, the same as the m~ 1 (wg), where the we are the extrema of v(w). 
To make the computation invariant it is sufficient to search for the extrema of the hyperbolic normalized 
densities (1 — \z\^) 2 ^{z) and (1 — |w| 2 ) 2 ^(u>) (which are invariant to Mobius change of coordinates). 
In our algorithm we consider the collection of Mobius transformations m = m(z; 9,p, q) defined by m{p) = q 
for every pair (p, q) € Is x Is' i an d every angle 9 € [0, 2tt). In order to compute the Mobius transformations 
in practice between two surfaces, we use the algorithm described in [111 112] . Furthermore, we discretize 9: 
9 = 2ir k/K, k = 0, 1, 2, K — 1. From every candidate Mobius m(z; 9,p, q) we build a candidate correspon- 
dence map C : S — > S' by the steps described in the next two subsections, deforming it to an area-preserving 

C m = f m O 771. 

One additional remark is that in the above algorithm we also consider all possible anti-Mobius transforma- 
tions m = m(z; 9,p, q) by taking fh(z) — m(z), where m is a Mobius transformation, and such that fh(p) = q 
for every pair (p, q) £ Is X Is>, and every angle 9 £ [0, 2n). 



4.2. Projection onto A(S,S'). Our goal now is to construct a map f m : <S — > S' such that f m o m E 
A(S,S'), and f m is (in some sense) "as close as possible" to the identity. 

Denote, as before, by n(z),v(w) the densities of the surfaces S, S' (resp.) over the unit disk V, defined by 
((mo$),j)[z] = [i(z)dzdz, and [w] — v(w) dw dw. Then, a simple and natural approach to define 

f m is via a "linear interpolation of the measures" technique due to Moser [T3] . The key idea is to look at the 
linear interpolant Ct : t > (1 — t)fi + 1 v, t G [0, 1] and to find a corresponding family of diffeomorphisms $ 4 
such that ($i)*d/x = dq t . (Here, as before, the (■)* notation, applied to a measure, means "push-forward", 
i.e. f*dfj, = dv is equivalent to the requirement that, for every Borel set f2, i/(/(J7)) = /u(fi).) Then the 
projection is defined as f m :— 4> m = $x- 

Dacaronga and Moser [2] used this strategy to construct an area-preserving map that takes a given density 
/ to a constant density. We will slightly generalize their formulation to achieve an area-preserving mapping 
<f> m :T> — > T> taking the area element dfj, to dv, that is 

(4.1) (0m)*e?M = dv. 

Other researchers have used Moser's technique to construct an initial guess in the further elaboration of 
an area-preserving map that would be optimal in the sense of mass-transportation cost [5]. Although 
Monge's mass-transportation provides a elegant way to construct correspondence maps, we believe that, 
because Euclidean (or hyperbolic) distances in the uniformization plane have no intrinsic meaning for the 
geometry of the problem, using them in the present context will not give a more meaningful answer than the 
straightforward result of Moser's procedure. More meaningful would be to use the surfaces' induced geodesic 
distances in a mass-transportation approach, but this is a much more challenging project, which we intend 
to tackle in future work. 

Since our measures are absolutely continuous w.r.t to the Lebesgue measures dz, dw respectively, we write 
dji = fj,(z)dz, dv — v(w)dw, using the conformal factors fi(z),v{w) as densities. Using the standard change 



of variables formula we see that (4.1) can be rewritten, in terms of the densities, as 



(4.2) " (<M*)) det(V0 m ) = m(z). 



We will be interested in a solution to (4.2 1 that is a diffeomorphism (f> m of V onto itself; in particular points 
on the boundary of the unit disk should be mapped to the boundary again. For the remainder of this 
subsection we will drop the subscript on cf> m , writing it as <fi for brevity. 
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Adapting Dacorogna and Moser's procedure [5] we define the diffeomorphism <f> by integrating, for t g [0, 1], 
a special time dependent vector field v t (z) (to be defined below): 

(4-3) ^U t (*)=M$ t (z)), 



dt 



(4.4) 



*o(*) 



for alii > 0, z E V 
for all z E T>. 



The desired map <p is then the end result of the integration, <j>(z) = $i(z). The vector field Vt is defined in 
three steps, as follows. We start by solving a Poisson equation with Neumann boundary conditions, 



(4.5) 
(4.6) 



0. 



in V 
on dT> . 



Aa = (j, — v, 
da 
dn 

[Note that, unlike Dacorogna and Moser we do not require <j>{z) = z for z € dT>; we impose only that the 
boundary of T> be mapped to the boundary - hence the use of Neumann instead of Dirichlet boundary 
conditions.] Next, a time- independent vector field v is defined by setting v(z) = Va(z). In the third step, 
we define the time- dependent vector field vt as 

v(z) 



vt(z) 



t ■ v{z) + (1 - t) ■ n(z) ' 



Establishing that <p{z) = $>x(z) provides a solution to (4.2) can be done by adapting Dacorogna and Moser's 
original proof. For completeness let us briefly describe the argument. First, we define an auxiliary function: 

(4.7) \{t,z)= (detV* t (z)) (t.i/($ t (z)) + (l-t)-A«(*t(*))); 
as we shall see below, this function satisfies 

(4.8) JU(M)=0. 



For the time derivative of the first factor we refer to [2]: 



(4.9) 



^ (det V$ t (z) ) = det V* t (z) ■ div v t ($ t (z)) . 



Differentiating (4.7) w.r.t. time gives thus 
(4.10) J^A(t, s) =det V* t • div v t ($ t ) • (t ■ v{$t 

+ det V® t (y( 



(l-t)- M ($ t )) 



/i(* t ) + (tvi/(* t ) + (i - t)VA*(*t), — $ t ; 



By the definition of w t we obtain 

div v = (div v t ) (t ■ v + (1 - t) ■ fi) + (<V^ + (1 - t)V(i, v t ). 
Together with (4.3) this leads to several cancellations in (4.10), resulting in 

d 



(4.11) 



dt 



\{t,z) = detV$ t (div v(<S> t 



Since v is defined as Va, and a satisfies (4.5), this implies (4. 



K$ t )-/i($ t ))). 

Therefore, 



A(0,z) = A(l,z). 

Because $o( z ) = z i we have A(0, z) = [J,(z), so that we have shown that 

fj,(z) = det V$i(z) v($i(z)). 



Finally, it is clear from the Neumann boundary conditions (4.6) that the vector field v(z) and therefore v t (z) 
is tangent to the unit circle at the boundary of the unit disk, that is (vt(z), z) = for all z 6 dT> and t > 0. 
This property ensures that integral curves &t(z) for z E dT> will stay on the boundary of the disk for all 
times t > 0. 
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Implementation details: We used the Matlab™ pde toolbox for all steps. For the first step (solving 
the Poisson equation) we used a triangular mesh with regular mesh size. The two densities are taken to 
be piecewise constant on the elements, with constants given by evaluating fj, and v at the midpoints of 
the triangles, providing the right-hand side of the PDE. Since the solution a of the PDE is also piecewise 
constant on the mesh elements, its gradient held v = grad a can be determined on each node of the mesh. By 
a nearest neighbor interpolation we approximate v as piecewise constant on the elements and use this to solve 
the ODE in the second step. This is done with a 4-stage Runge-Kutta method. Another implementation 
detail is that we add a small constant to the densities to avoid numerical inabilities for densities that have 
a minimal value close to zero. 

4.3. Thin-Plate Splines deformation. From a practical point of view we found it desirable to define our 
projection map as a composition of two maps: f m = <p m o £ m , combining the Moser map <p m defined above 
with a a preliminary smooth planar deformation £ m . The map Cm is used to locally align the peaks and 
valleys, already brought close together by the Mobius transformation m. Since we assume the two surfaces 
have equal (unit) area, i.e., J s dvols(x) = 1 = J s , dvo\s'(y), improving the alignment of peaks and valleys 
of the densities leads to less area distortion. This quick-and-dirty approximation jumpstarts the transition 
towards an exact area-preservation; although true area-preservation is achieved only after the second step of 
the deformation, an initial alignment by means of Q m removes some of the "workload" on <j) m . 

For the smooth deformation ( m , we use Thin-Plate Splines (TPS). In a first step, we label the points in 
Is and Is< as follows. We first apply m to the set Is, determine mutually closest points (with respect to 
the hyperbolic distance function) for the two sets m(Is) and Is>, and label them correspondingly, so that 
(PjiQj) £ nT- (Is) x Is 1 ,] = 1, denote the mutually closest pairs. In other words, we have 

d H (pj ,qj) < mm < min d H (pj,q), min d H (p, q 3 ) 

where the hyperbolic distance is dn{p, q) = tanh 

Next, we carry out a change of coordinates that maps the unit disk to the whole plane, by setting x( z ) — 
atan(|z|) zj \z\ with the inverse x -1 ^) = tan(|z|) z/ \z\. Set Pj = x{Pj)i Qj = xiVj), 3 = 1, We 
construct a thin-plate spline function £ m interpolating the Pj and Qj in the complex plane, i.e., Cm(Pj) = Qj- 
More explicitly, 

Cm [z] = X -1 ° TPS m o x, 

where 

n 

TPS m (z) — a Q + aiz + a 2 z + ^ 6 l T(|z - Pj\), 

i=l 

and T(r) = r 2 log(r). The coefficients cij,bi, j — 0, 1, 2,i — l,..,n are computed in the standard way by 
solving an (n + 3) x (n + 3) linear system [18] that imposes TPS m (Pj) = Qj, j = 1, . . . , n. "Sandwiching" 
TPS rn by the coordinate transformation \ guarantees that £ m takes the disk T> onto itself. 

4.4. Numerical experiments. Figures [3] and [4] demonstrate different aspects of the behavior of the algo- 
rithm described in the earlier sections. 

In the companion paper [4] an extensive analysis is performed for three biological data-sets, comparing the 
results of several algorithms to define the (dis)similarity between surfaces with those obtained by human 
experts. One of the methods illustrated in [3] uses the algorithm described here, and we refer the interested 
reader to that paper for many more figures and results. (In the interest of full disclosure, we confess that 
in many of the examples in 0] that used continuous Procrustes distances, we skipped the last step in Pr: 
the combination of an optimal Mobius transformation and TPS already gave results that were very close to 
area-preserving, and sufficed for the application at hand, so that we could skip the more time-consuming 
Moser transformation.) 
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Figure 3. Mesh deformation under Moser's procedure. Left: a "jiggled" initial 
mesh; Right: its deformation after solving the PDE on this mesh. 




Figure 4. The different components of the map. The left column shows the surface 
<S' (top, with colored squares enabling the viewer to track where different portions of the 
surface S are mapped), and the corresponding conformal factor v on the disk T>. The 
different steps of the algorithm "at work" in constructing the map C from S to S' are shown 
to the right of the vertical black line. From left to right: the optimal Mobius transformation 
to, which gives a uniformization of 5, with conformal factor \i\ the optimal alignment of the 
peaks in fi, via TPS, with those of v\ the transformation into a truly area-preserving map 
via Moser's technique. 
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